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Abstract 

The worldwide decline in honeybee colonies during the past 50 years has often been linked to the spread 
of the parasitic mite Varroa destructor and its interaction with certain honeybee viruses carried by Varroa 
mites. In this article, we propose a honeybee-mite-virus model that incorporates (1) parasitic interactions 
between honeybees and the Varroa mites; (2) five virus transmission terms between honeybees and mites at 
different stages of Varroa mites: from honeybees to honeybees, from adult honeybees to phoretic mites, from 
honeybee brood to reproductive mites, from reproductive mites to honeybee brood, and from honeybees to 
phoretic mites; and (3) Allee effects in the honeybee population generated by its internal organization such as 
division of labor. We provide completed local and global analysis for the full system and its subsystems. Our 
analytical and numerical results allow us have a better understanding of the synergistic effects of parasitism 
and virus infections on honeybee population dynamics and its persistence. Interesting findings from our work 
include: (a) Due to Allee effects experienced by the honeybee population, initial conditions are essential for 
the survival of the colony, (b) Low adult honeybee to brood ratios have destabilizing effects on the system, 
generate fluctuated dynamics, and potentially lead to a catastrophic event where both honeybees and mites 
suddenly become extinct. This catastrophic event could be potentially linked to Colony Collapse Disorder 
(CCD) of honeybee colonies, (c) Virus infections may have stabilizing effects on the system, and could make 
disease more persistent in the presence of parasitic mites. Our model illustrates how the synergy between the 
parasitic mites and virus infections consequently generates rich dynamics including multiple attractors where 
all species can coexist or go extinct depending on initial conditions. Our findings may provide important 
insights on honeybee diseases and parasites and how to best control them. 
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1. Introduction 

Honeybees are the world’s most important pollinators of food crops. It is estimated that one third of food 
that we consume each day mainly relies on pollination by bees. For example, in the United States, honeybees 
are major pollinators of alfalfa, apples, broccoli, carrots and many other crops, and hence are of economic 
importances. Honeybees have an estimated monetary value between $15 and $20 billion dollars annually as 
commercial pollinators in the U.S [22]. There are growing concerns both locally and globally that despite 
a 50% growth in honeybee stocks, the supply cannot keep up with the over 300% increase in agricultural 
demands [63] . Therefore, the recent sharp declines in honeybee populations have been considered as a global 
crisis. The most recent data from the 2012-2013 winter has shown an average loss of 44.8% of hives in the 
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U.S., and a total of 30.6% loss of commercial hives [53]. Some beekeepers have reported a loss of 90% of 
their hives [18, 37]. 

Between 1972 and 2006, the wild honeybee populations declined severely and are now considered vir¬ 
tually nonexistent [36, 61]. Hence the use of commercial honeybees for pollination is extremely important. 
Beginning in 2006, beekeepers began to report an unusual phenomenon in dying bee colonies. Worker bees 
would leave the colony to forage and never return, leaving the queen and the young behind to die. No 
dead worker bees were found at the nest sites; they simply disappear [12, 56]. This phenomenon is known 
as Colony Collapse Disorder (CCD), which is a serious problem threatening the health of honey bees and 
therefore the economic stability of commercial beekeeping and pollination operations. 

The exact causes and triggering factors for CCD have not been completely understood yet. Researchers 
have proposed several possible causes of CCD including stress on nutritional diet, harsh winter conditions, 
lack of genetic diversity, exposure to certain pesticides, diseases, and parasitic mites Varroa destructor which 
are also vectors of viral diseases of honeybees [22, 42], Even before CCD was detected in honeybee colonies, 
studies showed that most of the loses could be generally attributed to two main causes: the vampire mite, 
Varroa destructor, which feeds on host haemolymph, weakens host immunity and exposes the bees to a 
variety of viruses, and the tracheal mite, which infests the breathing tubes of the bee, punctures the tracheal 
wall and sucks the bee’s blood and also exposes the bee to a variety of viruses [45, 51, 32]. Since then, 
Varroa mites have been implicated as the main culprit in dying colonies. For example, in Canada, Varroa 
mites have been found to be the main reason behind wintering losses of bee colonies [20] , and more generally 
studies have shown that if the mite population is not properly controlled, the honeybee colony will die [50]. 
Recent studies also suggest that the Varroa mite could be a contributing cause of CCD since they not only 
ectoparasitically feed on bees, but also vertically transmit a number of deadly viruses to the bees [28, 27]. 
There have been at least 14 viruses found in honeybee colonies [4, 28], which can differ in intensity of impact, 
virulence, etc. for their host. For example, the Acute Bee Paralysis Virus (ABPV) affects the larvae and 
pupae which fail to metamorphose to adult stage, while in contrast the Deformed Wing Virus (DWV) affects 
larvae and pupae, which can still survive to the adult stage [57]. 

Mathematical models are powerful tools that could help us obtain insights on potential ecological processes 
that link to CCD, and important factors that contribute to the mortality of honeybees. Few sophisticated 
mathematical models of honeybee populations have been previously developed. DeGrandi-Hoffman et al. 
[14] produced the first time-based honeybee colony growth model. Martin [31] developed a simulation model 
consisting of ten components, which linked together various aspects of mite biology using computer software 
(ModelMaker); and Martin [32] later extended this model by including a bee model adapted from [14] to 
explain the link between the Varroa mite and collapse of the host bee colony. Wilkinson and Smith [62] 
proposed a difference equation model of Varroa mites reproducing in a honeybee colony. Their study focused 
on parameter estimations and sensitivity analysis. Simulation models are useful but may be too complex to 
study mathematically and obtain general predictions. 

More recently, mathematical models have been formulated to explore potential mechanisms causing CCD 
to the honeybee colony, it would be helpful to include the VARROAPOP model here because you can later 
compare the behavior and trends from your model with previously published models. Sumpter and Martin 
[54] modeled the effects of a constant population of Varroa mites on the brood and on adult worker bees, and 
found that sufficiently large mite infestations may make hives vulnerable to collapse from viral epidemics. 
Eberl et al. [16] developed a model connecting Varroa mites to CCD by including brood maintenance terms 
which reflect that a certain number of worker bees is always required to care for the brood in order for them 
to survive. They found an important threshold for the number of hive worker bees needed to maintain and 
take care of the brood. Khoury et al. [25, 24] developed differential equations models to study different death 
rates of foragers and the impact it had on colony growth and development. They then linked their results to 
CCD. Betti et al. [7] studied a model that combines the dynamics of the spread of disease within a bee colony 
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with the underlying demographic dynamics of the colony to determine the ultimate fate of the colony under 
different scenarios. Their results suggest that the age of recruitment of hive bees to foraging duties is a good 
early marker for the survival or collapse of a honeybee colony in the face of infection. Kribs-Zaleta et al. 
[26] created a model to account for both healthy hive dynamics and hive extinction due to CCD, modeling 
CCD via a transmissible infection brought to the hive by foragers. Perry et al. [41] examined the social 
dynamics underlying the dramatic colony failure with an aid of a honeybee population model. Their model 
includes bee foraging performance varying with age, and displays dynamics of colony population collapse 
that are similar to field reports of CCD. These models, no doubt, are insightful and provide us a better 
understanding on the potential mechanisms that link to CCD. However, most of these models only account 
for the honeybee population dynamics with mites or viruses but not both. 

The host-parasite relationship between honeybees and Varroa mites has been complicated by the mite’s 
close association with a wide range of honeybee viral pathogens. In order to understand how Varroa mite in¬ 
festations and the related viruses transmitted to honeybees affect honeybee population dynamics, and which 
may link to CCD, there is a need to develop realistic and mathematically tractable models that include 
both mite and pathogen population dynamics. The goal of our work is to develop a useful honeybee-mite- 
virus system to obtain better understanding on the synergistic effects of honeybee-mite interactions and 
honeybee-virus interactions on the honeybee populations dynamics, thus develop good practices to control 
these parasites to maintain or increase honeybee population. The most relevant modeling papers for our 
study purposes are by Sumpter and Martin [54] , and Ratti et al. [42] whose work examined the transmission 
of viruses via Varroa mites, using the susceptible-infectious (SI) disease modeling framework with mites as 
vectors for transmission. However, Sumpter and Martin assumed that the mites’ population is constant 
while Ratti et al. took no account of the fact that virus transmissions occur at different biological stages of 
Varroa mites and honeybees. 

In this article, we follow both approaches of Sumpter and Martin [54] and Ratti et al. [42], and propose a 
honeybee-mite-virus model that incorporates (1) parasitic interactions between honeybee and Varroa mites; 
(2) different virus transmission terms that account for the virus transmission among honeybees, between 
honeybees and mites at different stages of Varroa mites; and (3) Allee effects in the honeybee population 
generated by the internal organization of honeybees, including division of labor. Our proposed model will 
allow us explore the following questions: 


1. What are the dynamics of a system only consisting of honeybees and the disease? 

2. What are the dynamics of a system only consisting of honeybees and Varroa mites? 

3. What are the synergistic effects of Varroa mites and the disease on the honeybee population, and how 
may these synergistic effects contribute to CCD? 

4. How can we maintain honeybee populations? 


The structure of the remainder of the article is organized as follows: In Section 2, we first provide 
the biological background of honeybees, Varroa mites, and the associated virus transmission routes in the 
honeybee-mite system; then we derive our Si-type model for honeybees co-infected with the mite and virus. 
In Section 3, we perform local and global analysis of the proposed model and the related subsystems. The 
results from the analysis are then connected to biological contexts and implications. Additionally, we also 
explore numerical simulations of the subsystems and the full system to obtain the effects of each parameter 
in our system. In Section 4, we summarize our results and the related biological implications of our studies 
in finding potential causes of Colony Collapse Disorder. We also provide potential projects for future work. 
The detailed mathematical proofs of our theoretical results are provided in the last section. 
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2. Biological background and model derivations 

Honeybee colony: During the spring and summer, a honeybee colony typically consists of a single 
reproductive queen, 20,000 - 60,000 adult worker bees, 10,000 - 30,000 individuals at the brood stage (egg, 
larvae and pupae) and up to hundreds of male drones. During the winter, the colony typically reduces in size 
and consists of a single queen and somewhere between 8,000 - 15,000 worker bees [32], A large population 
of workers carry out the tasks of the bee colony, which include foraging, pollination, honey production and, 
in particular, caring for the brood and rearing the next generation of bees. The queen is the only fertile 
individual of the colony and has an average life span of 2 - 3 years [54]. During the peak season (in the 
summer), the queen lays up to 2000 eggs per day, where fertilized eggs produce female worker bees, or much 
more rarely queens, while drones develop from non-fertilized eggs [7]. The bees go through the following 
stages in development: egg (about 3 days), larvae (about 7 days), pupae (about 14 days), and adult. The 
life span of an adult worker bee also depends on the season. Workers usually have a lifespan of 3 - 6 weeks 
during the spring and summer, and are reported to live as long as 4 months during the winter [38]. The 
adult drone life span is typically 20 - 40 days, with reports of drone living up to 59 days under optimal 
colony conditions [38, 21]. 


Let Nh(t) be the total number of honeybees in the colony, including the larvae, pupae and adult bee 
(both hives and foragers) at time t. The subscript h means honeybees for all future notations. Define 
£h € [0,1] as the percentage of adult honeybee population in the colony, then (1 — £h)Nh is the brood 
population and is the ratio of adult honeybees to the brood in colony. Empirical study shows that 

the successful honeybee colony should have > 2 (see [47]). We should expect that the value of Ut 
varies with time. In our current model, instead of employing the explicit age structure model, we let be a 
parameter. Though brood and adult numbers change throughout the year, the ratio remains fairly constant 
and is a limiting factor in proportion of eggs that are reared into larvae and emerge as adults. An addition 
justification of this simplification is that the ratio of adult bees £/,, is indeed a constant at the steady state of 
the explicit age structure model. Under this simplification, we are able to obtain analytical results on how 
£h affects dynamics of the model with essential biological components. In the absence of mites and virus, 
the population dynamics of Nh(t) is described by the following nonlinear equation: 


K 


r(q h N h y- 

I< + {Z h N h f 


dhNh 


(1) 


where ' is the sign of the derivative with respect time; the parameter r is the maximum birth rate, specified 
as the number of worker bees born per day; the parameter \/~K is the size of the bee colony at which the 
birth rate is half of the maximum possible rate; and dh is the average death rate of the worker honeybees. 
The term h yi describes that the successful survival of an egg which will develop into a worker bee 

needs the care of adult honeybees (£ hNh ) inside the colony and also needs food brought in by the honeybee 
foragers. This term also implicitly assumes that more population of adult honeybees inside the colony can 
increase the survival of an egg developing into a worker which is supported by the empirical study showed 
in [47]. Our modeling approach of (1) follows the modeling idea in [16] for honeybee diseases and in [23] for 
the population of leaf-cutter ants. This term implicitly includes the internal organization of the honeybee 
population, such as division of labor. Model (1) implies that honeybee population is able to persist when £/,. 
is above an threshold, i.e., £/,, > 2dh More detailed analysis and related results are provided in the next 
section. 


Varroa mites: Varroa mites were first reported in Kentucky in the Bluegrass region of the Common¬ 
wealth in 1991 in U.S. They have since spread to become a major pest of honeybees in many states of 
U.S. [6]. Varroa mites are external honeybee parasites that attack both adult honeybees and brood, with 
a distinct preference for drone brood [39]. They suck the blood from both the adults and the developing 
brood, weakening them and shortening the life span of the bees which they feed on. Emerging brood may be 
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born with deformed wings. Untreated infestations of Varroa mites can cause honeybee colonies to collapse 
[30]. 

The mites go through a series of stages: larva, protonymph, deutonymph and then adult. Adult females 
undergo two phases in their life cycle, the phoretic and reproductive phases. During the phoretic phase, 
female Varroa feed on adult bees and are passed from bee to bee as they pass one another in the colony. 
During the phoretic phase, the female Varroa mites live on adult bees and can usually be found between the 
abdominal segments of the bees. The mites puncture the soft tissue between the segments and feed on bee 
lremolymph, harming the host [44, 8]. Mite reproduction can occur only if honeybee brood is available. A 
female mite enters the brood cell about one day before capping and will be sealed in with the larva. After 
the capping of the cell, it lays a single male egg and several female eggs at 30-hour intervals [58], and the 
mite feeds and develops on the maturing bee larva. When the host bee leaves the cell, the mature female 
mites leave the cell. The male mite dies after mating with his sisters, and if immature female mites are 
present they die as they come out of the cell, as they cannot survive once outside the cell. The adult female 
mite begins searching for other bees or larvae to parasitize. 

The phoretic period of the mite appears to contribute to the mite’s reproductive ability, which may last 
4.5 to 11 days when brood is present in the hive; or as long as five to six months during the winter when little 
or no brood is present in the hive. Consequently, female mites living when brood is present in the colony 
have an average life expectancy of 27 days, yet in the absence of brood, they may live for many months. 
In the average temperate climate, mite populations can increase 12-fold in colonies which have brood half 
of the year and 800-fold in colonies which have brood year-round. This period usually begins in late winter 
when brood rearing resumes from a winter period when little or no brood is present. The period of mite 
population increase continues through the spring and summer and peaks in the fall when brood rearing is 
nearly done. This makes the mites very difficult to control, especially in warmer climates where colonies 
maintain brood year-round [19]. 

Let N m (t) be the number of adult female Varroa mites in the honeybee colony in the absence of a virus 
where the subscript m means Varroa mites for all future notations. Varroa mites feed on the haemolymph of 
brood and adult honeybees, and their reproduction depends on the availability of the brood and the popula¬ 
tion of the reproductive female Varroa mites. Similarly to the case of honeybees, we incorporate an implicit 
age structure model of Varroa mites by defining a parameter £ [0,1] as the percentage of Varroa mites 
at the phoretic stage. This implies that the honeybee colony has the population size of (1 — £ m )-/V m as the 
reproductive Varroa mites and the population size of N m as the phoretic Varroa mites. This simplification 
still allows us to investigate the parasitic interactions between Varroa mites and honeybees rigorously with 
essential biological components. 

We model the parasitic interactions between Varroa mites and honeybees by using the Holling Type 
I functional responses, i.e. d(l — £/i)-/V/j(l — = aNhN m where a is the parasitism rate; the term 

(1 — £h)Nh is the honeybee brood population; the term (1 — £ m )-/V m is the reproductive female Varroa mites 
population; and a = d(l — ^/ t )(l — £ m ). Therefore, in the presence of Varroa mites N m , the dynamics of the 
honeybee population can be described as: 

K = khSI? - d h N h - aN h N m = - (4 + aN m )N h 

which implies that the parasitism decreases the life span of the honeybee, i.e., the average life span has been 
reduced to + * N after parasitism. Here we do not assume that parasitism would lead to the death of the 
honeybees for sure. The Varroa mite population depends on the nutrient obtained from honeybees, thus, 
the dynamics of Varroa mites population can be described as: 

N' m = caN h N m - d m N m 
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where the parameter a measures the parasitic rate of Varroa mites; c is the conversion rate from nutrient 
consumption obtained from honeybees to sustenance for Varroa mites reproduction; and d m is the natural 
death rate of Varroa mites. Therefore, in the absence of virus, the population dynamics of Varroa mites and 
honeybees can be described by the following two nonlinear equations: 

K = K r i%S^ -d h N h -aN h N m 

( 2 ) 

V m = CaNhNm dmNm 

The modeling approach of Model (2) has applied the traditional host-parasite modeling framework includ¬ 
ing non-lethal parasites [1, 2]. Model (2) also implies that Varroa mites population N m goes extinct if the 
population of honeybees Nh goes extinct. 


Varroa mites act as a disease-vector for virus transmissions: Varroa mites not only feed on 
host haemolymph and weaken host immunity, but they also expose honeybee colonies to at least 14 differ¬ 
ent viruses including deadly viruses such as ABPV and DWV. Varroa mites can transmit viruses in their 
reproduction phase to honeybee brood and during the phoretic phase to adult honeybees. To model the 
virus transmission between Varroa mites and honeybees during these two phases, we let Sh(t), S m (t) be 
the susceptible population of honeybees and Varroa mites, respectively; and be the virus in¬ 

fected population of honeybees and Varroa mites, respectively. Then the total population of honeybees is 
Nh(t) = Sh(t) + and the total population of Varroa mites is N m (t ) = S m (t) + 


The virus transmission between female Varroa mites and honeybees can occur in the following two phases 
of the Varroa mite life cycle: 

1. The honeybee colony has ^h.Sh susceptible adult honeybees; £,hlh virus infected adult honeybees, £, m Sm 
susceptible phoretic female Varroa mites; and £ m / m virus infected phoretic female Varroa mites. In 
the phoretic phase, female Varroa mites move between adult bees both spontaneously and just prior 
to the death of their host bee [54]. Following the approach of [32], we assume that virus transmission 
is frequency dependent, i.e., 


• We model the rate at which susceptible adult honeybees are virus infected by the infected phoretic 
female Varroa mites (IPFM) based on the approach of [32, 54, 16, 42]. This rate can be described 
as follows: 



x 


probability being infected after contacts 



X 


population of infected IPFM 


jhSh 

tihSh + £hlh 

probability contacts susceptible adult honeybees 


_ PmhtmShlm 

Sh+Ih 

which also implies that susceptible honey bees become virus infected at a rate proportional to the 
ratio of the population of the infected phoretic mites to the population of the total honeybees. 

• The rate at which susceptible phoretic female Varroa mites (SPFM) are virus infected by the 
virus infected adult honeybees (IAH) is: 



x 


probability being infected after contacts 



population of SPFM 


^h^-h 

tjhSh + tjhlh, 

probability being contacted by IAH 


_ PhmZm S m Ih 

Sh+Ih 
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2. The honeybee colony has (1 — £h)Sh susceptible honeybee brood; (1 — £h)Ih virus infected honeybee 
brood, (1 — £rn)Sm susceptible reproductive female Varroa mites; and (1 — £ m )/ m virus infected re¬ 
productive female Varroa mites. We do not assume that honeybee brood will die from parasitism, 
thus, parasitised honeybee brood or newborn Varroa mites will face virus infection if either honeybee 
brood or the reproductive female Varroa mite is virus infected. Chen et al. found that there is a direct 
relationship between virus frequency and the number of mites to which honeybee brood were exposed, 
i.e., the more donor mites that were introduced per cell, the greater the incidence of virus that was 
detected in the honeybee brood [10, 11]. This implies that the virus transmission rate between Varroa 
mites and the honeybee brood during the reproductive phase of mites is density dependent, i.e., similar 
to the term that describes the parasitic interaction between mites and honeybee. Therefore, we have 
as follows: 

• A newborn honeybee becomes virus-infected if it is parasitized by the infected reproductive female 
Varroa mites. Thus, the rate at which susceptible honeybee brood is virus-infected by the virus 
infected reproductive female Varroa mites (IRFM) is: 

fimh2 X (1 X d(l 

V«—^ ’ V ■ S - ’ V ■ 

probability being infected after contacts population of healthy honeybee brood parasitism by IRFM 


— (Imh2^ShIm- 

• The reproduction of Varroa mites depends on honeybee brood. The newborn Varroa mites become 
virus infected if either the brood or the female Varroa mites is virus infected. Thus, based on the 
formulation of the host-parasite interaction model (2), the rate at which infected newborn female 
Varroa mites (INFM) become virus infected depends on the parasitic interaction between mites 
and honeybees which can be described as ca [Ih(S m + Im ) + Shim] where the term Ih(S m + I m ) 
is the newborn Varroa mites infected with virus through virus infected honeybee brood; and the 
term Shl m is the newborn Varroa mites infected with virus through virus infected reproductive 
Varroa mites. 


The virus transmission among honeybees: The proportion of honeybees which can infect themselves 
is also dependent on the total number of susceptible and virus infected bees present in the colony, and hence 
frequency-dependent transmission is used [32], which is described as follows: 



x 


probability being infected after contacts 

h 

Sh + Ih 



the healthy honeybee population 

_ 0hS h I h 
Sh~\~Ih 


probability of contacting or being contacted by infected honeybees 


The reduced fitness of honeybees due to virus infections: The parasitic Varroa mites have been 
shown to act as a vector for a number of viruses including DWV, ABPV, Chronic Bee Paralysis Virus 
(CBPV), Slow Bee Paralysis Virus (SPV), Black Queen Cell Virus (BQCV), Kashmir Bee Virus (KBV), 
Cloudy Wing Virus (CWV), and Sacbrood Virus (SBV) [31, 32, 33, 27]. These virus infections contribute 
to morphological deformities of honeybees such as small body size, shortened abdomen and deformed wings, 
which reduce vigor and longevity, and they can also influence flight duration and the homing ability of 
foragers [27]. In our model, we assume that the virus infected adult honeybee population £/,//,. contributes 
to the reproduction of honeybees with a reduced rate p £ (0,1), therefore, the healthy honeybee population 
Sh can be modeled as follows, 
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PhShlh 

Sh+I h 


aSh{Sm + Im 


S'u = 


re h (s h + pi h f 

K + Ch (Sh + plh) ^°h ~t~ J-h^ 

reproduction of honeybees honey bee infected by themselves 

fimh (£ra Sh ) Im 


parasitism by mites 


• ( 3 ) 


S h + I h 


adult honeybees infected by the phoretic mites 


Pmh20tShIm. — dhSh 

honeybee brood infected by the reproductive mites 


And the virus infected honeybee population can be modeled by the following equation, 


Ih — Sh |] Sh + 'lb — OtIh(S m + Im) 

-( 


Consumed by mites 


(4) 


+ Ph )Ih 

natural honeybee mortality rate additional death due to virus infections 


Let /i m be the additional death rate of Varroa mites due to virus infections. Then the population of 
healthy Varroa mites S m and the virus infected Varroa mites I m can be described by the following set of 
nonlinear equations: 


S' = S n 


caSh — 


Ph m^m i h 

v ^ h I h 

the phoretic mites infected by adult honeybees 

Phm I h(im Sm) 


natural mortality rate of mites 


(5) 


I'm = ca [Ih(Sm + Im) + S h Im\ + ~ (<*™ + )Irr 

mites born with virus infections additional death due to virus infections 


Based on the discussions above, the full model of honeybee-mites-virus population dynamics is therefore 
modeled by the following system of differential equations: 


S' h = 

I'h = 
S'm = 
I' 


( _ . PhS h Ih _ l 3 mh (tmS h )Im 

K+e h (S h +pI h + ahDh S h +I h 8 h + I h 

Pmh2alm aSh{Sm V Im ) 

Sh |] Sh + Ih h" Pmh 20 lm^ ~ aIh(Sm + Im) ~ (dh + ph)Ih . 


( 6 ) 


Sm \caSh - l3h s ™‘^* h - dm] 

ca [Ih(Sm + Im) + Shim} + fihm g ~ {dm + Pm)Irr 


For convenience, let K = 4-, /3 mh = /3 m h£,m, Pmh = Pmh 2 Q., phm = Phm£,m- Then the full model (6) can 
























be rewritten as the following model 


S' h = 
I’h = 


where a > 0,c > 0, p £ [0,1] and the virus transmission rates f3h, $mh, fihm, Pmh £ (0,1)- In summary, 
the full honeybee-mite-virus model (7) incorporates (1) Allee effects of honeybees due to the cooperation of 
the internal organization; (2) parasitism interactions between honeybee and mites; (3) the vertical disease 
transmission mode modeled by the frequency-dependent disease transmission function during Varroa mites’ 
phoretic phase; (4) the horizontal disease transmission mode modeled by the density-dependent disease 
transmission function during Varroa mites’ reproductive phase; and (5) the reduced fitness of honeybees due 
to virus infections. The model (7) allows us to investigate the following scenarios: 

1. In the absence of the Varroa mites and virus, how population of honeybees may persist. 

2. In the absence of the virus, how the Varroa mites may affect the population dynamics of honeybees. 

3. In the absence of the Varroa mites, how virus infections may affect the population dynamics of hon¬ 
eybees. This case can apply to the situations that honeybees are infected by virus through ecological 
processes such as foraging other than parasitism by Varroa mites. 

4. In the presence of both the Varroa mites and virus infections, which conditions can lead to the extinction 
of Varroa mites, virus infections, and honeybees; and which conditions can guarantee the persistence 
of honeybee population. 

The rest of this manuscript will focus on the dynamics of Model (7) and the related subsystems. 


r(Sh+p!h ) 2 _ J C _ Ph.Sh.J-h _ Pmhtihlm _ f) C T _ rv C, /'C I T \ 
K+(S h +pI h ) 2 ah ^ h S h +I h S h +I h Pmh^h^m OO h {O m + l m ) 


0hShIh _ PmhSh 1 1 


Sh 

s m 


Phlh I Pmhlm I 5 T 
S h +I h ^ S h +I h ^ Pmh^m 


caS h - 1 hmI r h - d„ 


OIh(S m T ) (dh T dh)Ih 


( 7 ) 


s h +i h u ' m 


ccy. \ Ih{Sm + Im) + Shim ] + ^ +1™ ~~ Mm)-!) 


3. Mathematical analysis 

Let ^pj^| S(i=Jh=0 = 0 and s^+ Im \s m = Im =o = Define A' = {{S h ,I h ,S m ,I m ) £ R% : S h + I h > 
0 and S m + Im > 0}, then X can be considered as the state space of our model (7). To continue the 
analysis, let us define Nh = Sh + Ih, N m = S m + Im and N = cNh + N m as the population of honeybees, 
the population of mites, and the sum of the population of honeybees and mites, respectively. In addition, 
we let d = min{d/j, d m }, and define N* as the upper bound of the sum of the population of honeybees and 
mites and N c as the corresponding threshold, where 


N c = c 


_ J-\/(5) 2 -4A _ „ W(g) 2 - 

2 , — L 


4 K 


We let N£,N_h be the upper bound, lower bound of the population of honeybees, respectively, and N^,N^ 
be the corresponding thresholds, where 


Nh = 


-4 K 


N*h = 


VGfr) 2 - 4 * 


m = 


d h +^+«N*-\J { 4+^.K* ) -4K/P 2 d h+ ^+ a Jv*+d {d h +^+ aN * ) ~ 4 K/P 2 


N*h = 
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And we let S £ be the lower bound of the population of susceptible honeybees, and be the corresponding 
threshold, where 




+ 0 mh +*HN *^ 


d h+Ph + ^f^+0mh+<*)<. N *-cM-P 


-4 K 


S* h = 


d h.+f>h.+ f>TT $i’* +0mh+<*W*-c£LV 


\ 


d h+Ph+ f<m £*'* +(Pmh+°‘)( N *-cO-V 


-4 K 


\+\J{^Y- ik iv 


Define f b (x,y) = 3! ~ r ' v Kx> 2 -— and fb{x,y) = —— v Kx, 2 -—, then we have 


"\/(i) 2 - 4 ^ 


df b {x,y) „ n df b (x,y) ^ df b (x,y) ^ df b (x,y) 
ox oy ox oy 


which imply the following inequalities 


N c < Nf < Nf h <NZ< N* < N*/c. 


Theorem 3.1 (Basic dynamical properties). Assume that all parameters are strictly positive and p £ 
[0,1]. The model (7) is positively invariant and bounded in the state space X, which is attracted to the 
following compact set 

c = {{Sh, Ihi S m , I m ) £ Rjj, : 0 < c(Sh + Ih ) + (S m + An) = cNh + N m < IV*} 

2 > 2\/~ki and time is large enough. Moreover, the following statements hold for Model (7): 

> dh, then the total population of honeybees Nh is bounded by N£, i.e., 

limsup-ZV/^i) < N£. 

t—too 

> dit+Vh+aN an( ^ N h (Q) > N c h hold, then the total population of honeybees N b is persistent, 

N_h < liminf Nh(t) < limsup Nh(t) < < limsup N{t)/c = N*/c. 

t—too t—> oo 

• If the inequalities r /— > max [dh + fih + + {(3 m h + a)(N* — cNt), d ' t+Alh + aJV } with Nh{ 0) > 

2 y/K f — h p ) 

Sh( 0) > Sfr hold, then Sh is persistent with the following properties: 

< liminf Sh(t) < liminf Nh(t ) < liminf N(t)/c < limsup N(t) /c < N*/c. 

£->-oo t—>oo t—too t—too 


provided that 

•If-' 


2V K 


If- 


2V K 

i.e., 




The extinction equilibrium Eq = (0, 0,0, 0) is always local stable. Moreover, the system (7) converges 
to Eq globally if dh > —W holds; and it converges to Eq locally if the initial population satisfies either 


2v K 


JV(0) < N c or N h ( 0) < Nf. 


Notes: The positive invariance and boundedness results from Theorem 3.1 imply that our model is well- 
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defined biologically. In addition, Theorem 3.1 indicate follows: 


1 . 

2 . 

3. 


Initial conditions are important for the persistence of honeybees. 


The inequality 


i\JH 


> dh is a necessary condition for honeybee persistence, i.e., the 


growth rate r, small half saturation K , and the small death rate of honeybees dh- 


The small values of disease transmission rates /3h, Pmh, Smh', and small values of mite 
a are also important for the persistence of the healthy honeybee population Sh- 


large intrinsic 


attacking rate 


Recall that d = min{c4, d m }, N c = c- 




4 K 


m = 


dh+hh+<*N* 


V(i) 2 ~ 


, N* = c * ™ and 


Theorem 3.1 implies that under proper initial conditions, honeybees can persist if — 7 ^= > dh +P h + aN and Nh(0) > 

2 v k p 

N_h holds. Notice that ^ is the ratio of adult honeybees in the colony, and one of sufficient conditions that 
guarantee the persistence of honeybee population is the following inequality: 


r _ r£ h d h + Ph + aN* ^ r£ h _ a \Jid) d h + p h + ajg 

2 yfk ~ 2 Vk P 2 VK 2 p p 

Thus, Theorem 3.1 provides a critical function of the hives population ^ such that honeybee population 

t “1/(3) _ 4 p‘ 

can persist. In addition, notice that --—^- - is an increasing function of £/i, this implies that the 

higher hives to brood ratio £h, the better honeybee growth will be, and more likely persistent for honeybees. 
This is supported by the empirical study of [47]. In the following theorem, we provide theoretical results 
on the sufficient conditions that lead to the persistence and the extinction of disease population or Varorra 
mites population. 


Theorem 3.2 (Persistence and extinction of disease or mites). The following statements hold 
• If N* < =**■, then the total population of mite N m goes extinct, i.e., 

limsup./V m (f) = 0 

t—t OO 


where system (7) is attracted to the mite-free invariant set MF = {(Sh, Ih, S m , Im) £ R.+ : S m + I m = 
0}, and its dynamics is equivalent to the following two-D model ( 8 ) 


S' h 


I'h 


_ r(Sh+pIh ) 2 _ J C _ C Phlh 

K+{S h +pIh ) 2 dhbh bh S h +I h 

= S h ^-(d h + p h )I h 


( 8 ) 


If 


2V K 


> dh, and Sh( 0) > S%, then the total population of honeybees persists while the 


healthy mite population S m goes extinct, i.e., 


lim sup S m ( t ) = 0 

t—> OO 
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where the system (7) is attracted to the healthy-mite-free invariant set HMF = {(Sh, Ih, S m , I m ) £ 
: S m = 0} and its dynamics is equivalent to the following three-D system (9): 


S'h 


I'h 


1~(Sh+pIh) 2 J C PhShlh fimhSh 1 1 

K+(S h+P I h ) 2 ahDh S h +I h S h +I h 


IL = 


s h 

cal, 


Phlh I Pmhlrr 


Sh~\~Ih Sh~\~Ih 

h 


fimhShlm & Sherri 

(dh ~h fJ'h'jlh 


(9) 


c dm+Ll-n 

”77 


Assume that — T -j= > dh + Ph + —b (f3 m h + a)(N* — Nff) and Sh( 0) > Sf. Then the disease 


2 yj K 

min < P h ,cP mh +cP mh +caS^ > 

I = clh + I m persists if the inequality -> 1 holds. 

max < ( d h +p h ),(d m +p m ) > 

Assume that —'-/= > dh +^ h + a x anc [ NhlO) > N_h■ Then the disease I = clh + I m goes extinct if 

K p 


the following inequality holds 


max s 0h+ ^ h ^f r ,cPmh+cPmh+caNZ > 

_L- =h - L <- y 


Under this condition, the system 


^(d h +p h ),(d m +p m )^ 

(7) is attracted to the virus-free invariant set DF = {(Sh, Ih, Sm, Im) £ '■ h + I m = 0} and its 

dynamics is equivalent to the following two-D model (10) 


S{ = 


rSt 


h K+S 2 


dhSh aShSr, 


( 10 ) 


S 7n — caShSm d m S„ 


Notes: The results of the reduced dynamics in Theorem 3.2 can be easily obtained by the theory of asymp¬ 
totically autonomous systems [9]. The detailed proof of our results are provided in the last section. 

^ + irry— 

According to Theorem 3.2, the condition N* = c— —-—^- < — can lead to the extinction of the mite 

population. Therefore, we can conclude that large values of the death rate of mites, d m , can lead to the 
extinction of the whole colony; and large values of the death rate of mites d m , small values of mite attacking 
rate, a, and its energy conversion rate c, can lead to either the extinction of the whole mite population N m 
or the extinction of the healthy mite population S m . Here we would like to point out that it is possible to 
have the persistence of virus infected mites while the healthy mite goes extinct (see the resulting dynamics 
(9) when the healthy mite goes extinct). In addition, the results of Theorem 3.2 also suggest that: 1. the 
persistence of the virus requires a large value for the disease transmission rate between adult honeybees, /3 h, 
or that the disease transmission rate between honeybee brood and reproductive mites, fd m h', or small values 
of total death rates of honeybees, dh + p>h, and mites, d m + /z m ; 2. the extinction of the virus requires small 
values of all disease transmission rates, i.e., small values of (3h, Pmh, Pmh, $hm', or large values of total death 
rates of honeybees and mites. 

Theorem 3.2 provides sufficient conditions that the full system (7) reduces to the virus-free subsystem (10), 
the mite-free subsystem ( 8 ), and the healtliy-mite-free subsystem (9). In the following three subsections, we 
explore the global dynamics of these subsystems. 
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3.1. Dynamics of the virus-free subsystem: only parasitism by mites 

Theorem 3.2 in previous section suggests that either small values of all virus transmission rates or large 
values of total death rates of honeybees and mites can lead to the extinction of the virus infected honeybees 
and mites, which gives the following virus-free dynamics (10): 


S' h 


= l S c2 - d h S h - aS h S„ 


S m ca S h S fn d rn S n 


The dynamics of the virus-free dynamics (10) (i.e., the dynamics of the parasitism interactions between 
honeybees and Varroa mites) can be summarized by the following theorem: 


Theorem 3.3 (Dynamics of the virus-free subsystem). Let H* = — ( 
The virus-free subsystem (10) can have one, three, or four equilibria. 


M* = 


rH * 


- d h 


« Ik+(H*) 2 

The existence and stability 
tions for these equilibria are listed in Table 1. The global dynamics of the virus-free subsystem (10) 


condi- 
can be 


Equilibria 

Existence Condition 

Stability Condition 

(0,0) 

Always exists 

Always locally stable 

Wi,0) 

—jt > d h 

2v K 

Saddle if N£ < &* = H*; Source if Nf > ^ = H* 

W.0) 

—7w > d h 

2v K 

Sink if Nf < duL = H*-, Saddle if Nt > ^ = H* 

ri olc ’ h, olc 

(. H*,M*) 

m < ^ = h* < n* 

_ Lk __ lL- 

Sink if H* > Vl<; Source if H* < \Jk. 


Table 1: The existence and stability of equilibrium for the virus-free subsystem (10), where 

^r + \/(^r ) 2 ~ 4 


and H* = 4bl, M* = 1- \ - rH " , 

olc ’ a [k+(U *) 2 


-d, 



N h = 


summarized as follows: 

1. The system (10) converges to extinction (0,0) for almost all initial conditions if — 'j= < dh or < Nf. 


2 vk 


2. If N£ < , depending on initial condition, the trajectory of (10) converges to either (0,0) or (iV^,0). 

3. If Nf < dff < N* then Model (10) has a unique interior equilibrium ( H*,M *) which is locally asymp¬ 
totically stable when — > yffc and is a source when — < yffc. 


Notes: Theorem 3.3 provides us a global picture on the dynamics of the virus-free subsystem (10), i.e., the 
honeybee colony only virus infected with mites but not the virus. By applying the results in [55, 59], we 

can conclude that the virus-free subsystem (10) undergoes a subcritical Hopf-bifurcation at ^ = \fH. The 

subsystem (10) has a unique unstable limit cycle around ( H*,M *) whenever — < \[k. In this case, the 
periodic orbits expand until it touches the stable manifold of the boundary equilibrium (Nf, 0) which leads 
to the extinction of both honeybees and the parasitic mites. We refer to this phenomena as a catastrophic 
event which could be linked to CCD. Our theoretical results also suggest that a small death rate for mites 
and a large parasitism rate can destabilize the system. 


Linking to CCD: To illustrate the catastrophic event, we use reasonable parameters from [54, 42] . Let the 
reproduction of egg per day during summer be r = 1500; and the population size of the honeybee colony 
at which the birth rate is half of the maximum possible rate be \[k = 2000; the natural death rate of 
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Figure 1: Time series (in days) of Model (10) when r = 1500, a = 0.005, c = 0.01, = 0.01, dm = 0.1: population of 

honeybees is in red while Varroa mites is in blue. 


honeybees is dh = 0.01; the parasitism rate is a = 0.005; the energy conversion rate is c = 0.01; and the 
natural death rate of mites is d m = 0.1. This set of parameter values gives ^ < \f~K which implies that 
a catastrophic event will occur (see Figure 1; the population of honeybees is in red and collapses around 
time =200 day). 

Stochastic effects and oscillations: Theorem 3.3 implies that if the inequality ^ < \f~K holds, 
then the virus-free subsystem ( 10 ) has a unique unstable limit cycle around ( H*,M*) where for all initial 
conditions either the system converges to ( 0 , 0 ) quickly or the system experiences the expanding oscillations 
leading to the extinction eventually. The oscillating extinction in the later type is driven by the deterministic 
dynamics. The extinction fate of the system cannot be prevented by introducing stochastic effects, however, 
introduced stochastic effects may cause the system going to extinction more quickly without expanding os¬ 
cillations. 

Note that K = ^ and a = d(l — 60(1 — Cm) where 6 i, Cm are ratio of adult bees and ratio of phoretic 

. h . 

stage of Varroa mites in honeybee colony, respectively. The catastrophic event occurs when 

d m ^ dm (1 — 60 dt( 1 — U)cVK 

- < VA w -< --- 7"-—r < - - -. 

ere d(l — Cm)cv K 61 (1 60 d m 

This inequality provides a critical low hive to brood ratio that can destabilize the system and cause the 
sudden extinction of honeybees. 


14 










3.2. Dynamics of the mite-free subsystem: only virus infections 

According to Theorem 3.2, if the honeybee population is too small, e.g., N* < —, then the dynamics of 
(7) is equivalent to the following mite-free dynamics (8) 


S' h 


= Ks-d h S h -S h J&, 


K+(S h +pI h ) 2 


Sh+Ih 


I'h 


= S h7 ^- h -(d h + f x h )I h 


To continue studying the dynamics of the mite-free dynamics (8) , we define a = —- as the ratio of 

d h+i*h ~ 1 

the susceptible honeybee population to the virus infected honeybee population; = d ^ff as the basic 
reproduction number, i.e., the number of secondary cases which one case would produce in a completely 
susceptible population; d = (a + 1 )dh + Hh = dh J + Th as the updated average death of the honeybee 

due to virus infections. In addition, we let = all,k = 1,2 and 


I\ = 




(p+p ) 2 


T 2 - 


(5 ) 2 




The dynamics of the nrite-free system (8) can be summarized by the following theorem: 

Theorem 3.4 (Dynamics of the mite-free subsystem). The mite-free subsystem (8) can have one, three, 
or five equilibria. The existence and stability conditions for these equilibria are listed in Table 2. In addition, 


Equilibria 

Existence Condition for Existence 

Stability Condition 

(0,0) 

Always exists 

Always locally stable 

W),0) 

- J T ^>d h 

2 \J K 

Saddle if 3 ?q <1; Source if Hq >1 

0) 

—7= > d h 

2v K 

Sink if 3^ <1; Saddle if > 1 

(Sill) 

> 1 and ^ > „t p 

Always a saddle. 

(Si ID 

> 1 ^ 2 Ifk " a +P 

Always locally asymptotically stable 


Table 2: The existence and stability of equilibrium for the mite-free subsystem (8). We have = —A—, a = 


(a + 1 )d h + Uh = d h 


+ V-h, N? - 


a7 , _ ^+\/fe) 2 - 4 ^ and ji = S~\/(S) 2 - 4 (T 


+p ) 2 


d = 


r 2 - 

1 h — 




SI, S jt = allk = 1,2. 


the global dynamics of the mite-free subsystem (8) can be summarized as follows: 

1. The trajectory of (8) converges to extinction (0,0) for all initial conditions in if one of the following 
conditions hold: 

= < d h , or 


2 \/K 


> 1 and d h <^< 


2V K 


2. The trajectory of (8) converges to either (0,0) or (IV)), 0) for almost all initial conditions in . 
inequalities IRq < 1 and —> dh- 


•2\J K 


3. The trajectory of (8) converges to either (0,0) or (S^, Iff) for almost all initial conditions in 


if the 
if the 


inequalities Hq > 1, 


> max 


2 \JK 




hold. 


15 












































Notes: Theorem 3.4 implies that the mite-free subsystem (8) has relatively simple dynamics, i.e., no limit 
cycle. The results show the following interesting findings: 


1. Honeybees can persist with proper initial conditions if the virus transmission rate among honeybees 
j3h is not large, i.e., IRq < 1. 

2. Both honeybees and the virus can coexist if (3h is in the medium range, i.e. 3 £q > 1 and < 


2V K 


3. However, the large virus transmission rate among honeybees {3h can drive honeybees to extinction. 
This occurs when the inequalities > 1 and dh < 


< _2_ hold. 

2 \! K a+p 


3.3. Dynamics of the healthy-mite-free subsystem 

According to Theorem 3.2, if ^ > dh, Nj) < jff, and 5^(0) > S r f : , then the total population of 
honeybees persists while the healthy mite population S m goes extinct, i.e., 

lim sup Sm(t) = 0 


Let 


S' h = 


r{Sh+pIh .) 2 _ J C _ QhShlh _ PmhShlm _ 

K+(S h +pI h ) 2 hDh S h +I h S h +I h 


where the system (7) is attracted to the healthy-mite-free invariant set HMF = {(Sh, Ih, S m , Im) £ 
S m = 0} and its dynamics is equivalent to the following three-D system (9): 

3 rnhSh^m G: S f h l rn 

Sh + Ih ^Sh+Ih fimhlm ~ <+IhI m ~ (dh + Ph)Ih 

I'm = Cal m \l h +S h 


*■+ ' 


I'h = S h 


r{^f^-i h + P i h y 


h(h) = 


f2(Ih) = 


«+( 


dm+t+m 


I h+p I h) 


'Z-dh dm fX m -Phlh 


d-m+t+rr 




a+ 


dm+Vn 




dm+F' 


‘t+hmh 


The dynamics of the healthy-mite-free subsystem (9) can be summarized by the following theorem: 


Theorem 3.5 (Dynamics of the healthy-mite-free subsystem). If — j— > dh and < 


dm + Pr, 


2V K 


■, then the 


population of virus infected mites goes extinct in the subsystem (9), i.e., 


lim sup I m ( t ) = 0 

£—>•00 


which reduces to the following mite-free model (8). In addition, the following statements hold 

1. If (Sh, Ih, Im) is an interior equilibrium of the healthy-mite-free subsystem (9), then Ih is a positive 

intercept of fi(h) and / 2 (4) subject to 0 < I h < tSiidJhs , ; S h = - Ih and I m = fi(h)- 

2. The healthy-mite-free subsystem (9) has no interior equilibrium if the inequality rca < dh(d m + Pm) 
holds. 


16 































3. Assume that pr > d h+Uh+aN an ^ iV/ l (0) > N %. TTien both virus infected honeybee population Ih 


and virus infected mites I m persist if the inequalities > 


dm + U-n 


and — 

u d h +fj. h 


< 1 hold. 


Notes: Theorem 3.5 implies that, under the condition of , the virus infected mite population 

I m goes extinct in the healthy-mite-free subsystem (9) which reduces to the mite-free subsystem (8) that 
we studied in the previous subsection. In addition, Theorem 3.5 shows that the subsystem (9) has no 
interior equilibrium if the inequality rca < dh(d m + p m ) holds. Therefore, we could expect the extinction 
of I m for small values of r,c,a and large values of dh, d m , p m . This has been confirmed by numerical 
simulations. The population of honeybees and virus infected mites in (9) experiences sudden collapse when 
we (1) increase the values of c, a, K, and the related virus transmission rates, or (2) decrease the values of 
dh, p,r, pL m . The biological implication for this dynamics is that increasing or decreasing the values of these 
parameters destabilizes the system and generates fluctuated dynamics. The destabilizing effects generate 
unstable oscillations. The amplitudes of oscillations increase until they touch the stable manifold of the 
extinction equilibrium, which cause the collapse of the whole colony. The destabilizing effects of c, a , K , d m 
can be explained through the dynamics of the virus free subsystem (10) that we have studied in Theorem 
3.3. 

• Decreasing the values of c, a, K can stabilize the system; small values of c, a can cause the extinction 
of the virus infected mite population / m , and lead to the coexistence of Sh and Ih- 

• Increasing the value of ph can stabilize the system but large values of ph can cause extinction of the 
whole colony due to the initial oscillations. 

• Decreasing p m can destabilize the system; while increasing it can stabilize the system; large values of 
p rn can lead to the extinction of I m and the persistence of Sh,Ih- 

• Decreasing the value of p could destabilize the system, thus causing the extinction of the colony. 

• Increasing the virus transmission rates (i.e., j3h , Phm, Pmh, fthm) can stabilize the system, while decreas¬ 
ing their values can destablize the system and cause the extinction of all species. 


Synergistic effects of parasitic mites and virus infections: If there is no mites in the system, ac¬ 
cording to Theorem 3.4, the mite-free subsystem (8) reduces to the only healthy honeybee population, i.e., 
the disease goes extinct whenever the initial population of honeybees is above Nf, p— > d h+Uh+aN , anc [ 

the basic reproduction number IRq = dl ^ /]ih < 1 (see two figures in the first row of Figure 2 where virus 
infected honeybees go extinct (the black curve in right) and the healthy honeybees persist (the red curve in 
red)). However, when the virus infected mites are in the colony, i.e., the healthy-mite-free subsystem (9), 
both disease and mites can persist under proper conditions (see Figure 3 where virus infected honeybees (the 
black curve), the healthy honeybees (the red curve), and the virus infected mites persist (the cyan curve) 
with the healthy mites going extinct (the blue curve). For example, the synergistic effects of parasitic mites 
and virus infections has been illustrated in Figure 2-3 when r = 1500; I\ = 1000000; p = 0.9; dh = .15; ph = 
0.1; a = 0.005; d m = 0.1; p m = 0.01; c = 0.005;/?h = .24; jd rn h = 0.03 ;j3 m h = .005 ; fihm = 0.03. These are 
reasonable parameter values derived from [54, 42]. 
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Figure 2: Population dynamics of the subsystems of the the honeybee-mite-virus model (7) when r = 1500; K = 1000000; p = 
0.9; dh = .15;= 0.1; a = 0.005; dm = 0.1; /x m = 0.01; c = 0.005; fth = .24;/3 m ^, = 0.03; fimh = .005;/3^ m = 0.03. The left 
figure in the first row is the healthy honeybee population (the red curve) and the right figure in the first row is the virus 
infected honeybee population (the black curve) in the mite-free subsystem (10) when Sh{ 0) = 4001, 7^(0) = 10. The left figure 
in the second row is the healthy honeybee population (the red curve) and the right figure in the second row is the healthy mite 
population (the blue curve) in the virus-free subsystem (8) whenS'^(O) = 4001,S' m (0) = 5. 


3-4- Dynamics of the full system 

Recall that the full system (7) of honeybee-mite-virus interactions can be described by the following set 
of equations: 


Q' 
D h 


r(S h +pI h ) 2 


dhSh 


PhSh Ih _ PmhShlrr 


K+(S h +pI h ) 2 ahDh Sh+Ih ' "s’k+ih" PmhShlm OtS h (S m + I m ) 


I' h 


O' 

J m 

r 


s h 

Sm 


Phlh | Pmhlm | a T 
S h +I h + S h +I h + Pmh^n 


C- _ Phmlh _ J 

cao h Sh+Ih a„ 


&Ih(Sm T Im) (dh T l^h)Ih 


— ca [Ih(S m + I m ) + Shim] + — (dm + Mm)X, 


The results from the previous section provide us a complete picture of the dynamics of the subsystems of 
the full system (7). In this subsection, we explore the dynamics of the full system as the following theorem. 


Theorem 3.6 (The persistence of honeybees). Assume that —> dh- If N* < ^ and Hq = < 

1, the full system (7) converges to the disease-mite-free set DMF = {(Sh, Ih, S m , Im) £ R+ : S m + Ih + Im = 
0} where the system (7) is reduced to the following one-D system (11).' 


St 


rSj 

k+si 


— dhSh 


whose dynamics can be summarized as follows: 


( 11 ) 
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Figure 3: Population dynamics of the honeybee-mite-virus model (7) when r = 1500; K = 1000000; p = 0.9; dh = .15; ^h = 
0.1; a = 0.005; dm = 0.1; = 0.01; c = 0.005; ^ = .24;/3 m/l = 0.03;= .005; P hrn = 0.03 and S h ( 0) = 4001,7^(0) = 

10, Sm{ 0) = 5,7 m (0) = 10. The healthy honeybee population Sh is in red; the virus infected honeybee population Ih is in black; 
the healthy mite population Sm is in blue; and the virus infected mite population 7 m is in cyan. 
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4- Assume that —> dh+un+aN an ^ N h IQ^ > fVS. Then the total mite population N m = S m + I m 
i\J k p 

persists if the inequalities and IRq = df+^. h < - 1 hold. 


Notes: Theorem 3.6 along with Theorem 3.1 - 3.3, we can conclude that the extinction of disease occurs 
when all values of all disease transmission rates, (3h, $mh, Pmh, Phm are small; with the consequence that the 
full system (7) converges to either (0,0, 0,0) or (iV^, 0, 0,0) when IRq 1 = ^ < 1 while (7) converges to either 
(0,0,0, 0) or (iJ*, 0, M*, 0) when 1 < The persistence of disease or mites indicates the persistence 

of honeybees even though the population may be low under the influence of disease or mites. Theorem 3.6 
provides a summary on sufficient conditions when honeybees can persist in the full (7) alone, with mites, 
or with disease. The item 4 of Theorem 3.6 is consistent with the results from Theorem 3.5 regarding the 
synergistic effects of parasitic mites and virus infections: If there is no mites in the system, according to 
Theorem 3.4, the mite-free subsystem (8) reduces to the only healthy honeybee population, however, in the 
presence of mites, both disease and mites in the honeybee-mite-virus system (7) can persist under proper 
conditions (see Figure 2-3 for more details). 

The dynamics of the full system (7) can be extremely complicated. We are unable to obtain an explicit 
form of the interior equilibrium and the related stability. We perform a series of numerical simulations to 
explore how different parameters affect the population dynamics. The effects of r, c, a, dh , d m , p m , /ah, p and 
virus transmission rates are similar to our observations for the subsystem (9). More specifically, we have the 
following observations: 

1. Effects of r: Increasing r can stabilize the system, but increasing it too large can drive healthy mites 
S m to extinction while the population of virus infected mites I m increases. Decreasing the values of 
r can destablize the system, and cause the extinction of mites. And too small values can cause the 
whole colony to become extinct. 

2. Effects of c: Increasing can destabilize the system and cause extinction of all species. Decreasing its 
value can stabilize the system but too small values can drive the extinction of mites. 

3. Effects of a: Increasing can destabilize the system and cause extinction of all species while decreasing 
can stabilize the system and increase the healthy honeybee population. Small values of a can drive 
the mite population to extinction, and too small values can drive both mites and virus to extinction, 
and only healthy honeybees are left. 

4. Effects of d m : Increasing d m can stabilize the system, and drive S m to extinction first. Increasing it 
further can lead to the extinction of mites, and the system approaches the limiting mite-free system 
(8). Decreasing can destabilize the system and drive the extinction of the virus. Too small values can 
cause the whole colony’s extinction. 

5. Effects of dh'- Increasing can stabilize the system. Large values can drive the virus extinctions, however, 
extremely large values can lead to the extinction of the whole colony. Decreasing can destabilize the 
system which may cause the extinction of the colony under certain conditions. 

6. Effects of p m - Increasing can drive virus extinction. Extremely large values can lead to the extinction 
of the whole colony. Decreasing can destabilize the system. Small values can drive the healthy mites 
extinct, and extremely small values may cause the extinction of the whole colony. 

7. Effects of ph' Increasing can cause the extinction of virus. Decreasing can stabilize the system, and 
small values can drive healthy mites to extinction. 

8. Effects of virus transmission rates: Increasing can destablize the system and cause the extinction of 
healthy mites S mi while extremely large values may drive all populations extinct. Decreasing can drive 
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the extinction of virus. 


3.5. Mechanisms of collapse dynamics and synergistic effects 

Let r = 1500; K = 1000000; p = 0.9; d h = .15; p h = 0.1; a = 0.005; d m = 0.1; p m = 0.01; c = 0.005; = 

.2A-,f3 m h = 0.03 ;/3 mh = .005 ',fhm = 0.03 (Figure 2-3) and r = 1500; p = 0.9; K = 1600001; dh = 0.15; ph = 
0.1; a = 0.05; c = 0.005; d m = 0.1 ;/x TO = 0.01; /3 h = 0.3; fd m h = 0.08; j3 mh = 0.001 ;j3hm = 0.03 (Figure 4-5). 
These are reasonable parameter values derived from [54, 42]. We use these two sets of parameters as illustra¬ 
tions to explore the synergistic effects of parasite mites and virus infections as well as potential mechanisms 
linking to CCD (see Figure 1 and Figure 4-5). These comparisons suggest the following: 


1. Synergistic effects of parasitic mites and virus infections: Based on the two sets of parameters, we have 
the following two typical scenarios: 

(a) Under the parameter values of r = 1500; IF = 1000000; p = 0.9 \dh = .15; ph = 0.1; a = 
0.005; dm = 0.1; p m = 0.01; c = 0.005; 0 h = .24; fi rnh = 0.03; (3 mh = .005; (3 hm = 0.03: 

If there is no mites, the mite-free system (8) (i.e., the dynamics of the healthy honeybee and 
the virus infected honeybee) converges to only the healthy honeybee with virus infected honey¬ 
bees go extinct (see the first row of Figure 2). 

If there is no virus, the virus-free system (10) (i.e., the parasitism dynamics of the healthy hon¬ 
eybee and the healthy mites) converges to a stable equilibrium where both the healthy honeybee 
and the healthy mites can persist (see the second row of Figure 2). 

However, if honeybees, mites, and virus are all presented in the system (i.e., the full system 
(7)), then both virus infected honeybees (black curve) and virus infected mites (cyan curve) can 
persist (see Figure 3). 

This implies that the presence of mite population can promote the persistence of disease. 

(b) Under the parameter values of r = 1500; p = 0.9; K = 1600001; dh = 0.15; ph = 0.1; a = 0.05; c = 
0.005; d m = 0.1; p m = 0.01; fi h = 0.3; /3 mh = 0.08; p mh = 0.001; p hrn = 0.03: 

If there is no mites, the mite-free system (8) (i.e., the dynamics of the healthy honeybee and 
the virus infected honeybee) converges to a stable equilibrium where both the healthy honeybee 
and the virus infected honeybees can persist (see the first row of Figure 4). 

If there is no virus, the virus-free system (10) (i.e., the parasitism dynamics of the healthy honey¬ 
bee and the healthy mites) converges to the extinction of both species through catastrophic event 
(see the second row of Figure 2). 

However, if honeybees, mites, and virus are all presented in the system (i.e., the full system 
(7)), then both honeybee and mites go extinct (see Figure 3). 

This implies that the presence of the unstable mite population can lead to the extinction of 
honeybees. 
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Figure 4: Population dynamics of the subsystems of the the honeybee-mite-virus model (7) when r = 1500; p = 0.9; K = 
1600001; dh = 0.15;= 0.1; a = 0.05; c = 0.005; dm = 0.1; p,m = 0.01; (3^ = 0.3;— 0.08; firnh = 0.001; fthm = 0.03. The 
left figure in the first row is the healthy honeybee population (the red curve) and the right figure in the first row is the virus 
infected honeybee population (the black curve) in the mite-free subsystem (10) when Sh{ 0) = 7684,7^,(0) = 1700. The left 
figure in the second row is the healthy honeybee population (the red curve) and the right figure in the second row is the healthy 
mite population (the blue curve) in the virus-free subsystem (8) when 5^(0) = 410, 5 m (0) = 35. 


2. Linking to CCD: In the absence of virus infections, the subsystem (10) goes through the catastrophic 
event which causes the extinction of honeybees (see the honeybee population in the red curve of Figure 
1). This property has been inherited by the full system (7) as honeybee population goes extinct sud¬ 
denly for most initial conditions (see the honeybee population in the red curve of the fist left figure in 5). 


Our analysis and simulations suggest that it is important to include both mites and virus in studying 
the population dynamics of honeybees as we proposed the full system (7) due to the synergistic effects of 
parasitism induced by mites and the virus infections as well as the catastrophic event from the parasitism 
interactions between mites and honeybees. In other words, if we only consider the honeybee versus virus 
dynamics as described by the subsystem (8), or only consider the parasitism interactions between mites and 
honeybees as described by the subsystem (10), we are not able to capture the full mechanics that can lead 
to the extinction of honeybees or the biological implications on the persistence of disease in honeybees. In 
addition, the full system (7) has rich dynamics which can possess multiple attractors. Thus, depending on 
initial conditions, the full system (7) can either experience extinction or have coexistence of both honeybees 
and mites. More sophisticated mathematical analysis is needed in order to understand the detailed dynamics. 


4. Discussion 

The association of virus infection with Varroa mites infestation in honeybee colonies causes great concern 
for researchers and beekeepers [11]. Many studies have suggested that Varroa mite infestations could be a 
key explanatory factor for the widespread increase in annual honeybee colony mortality and has been impli¬ 
cated as a contributing factor leading to CCD [35]. In this paper, we derive and study a honeybee-virus-mite 


22 



















500 


400 



10 

8 



300 

■C 

C/) 

200 


Honeybee-Mite-Virus Dynamics ( S h -I h -S m -I m ) 



Honeybee-Mite-Virus Dynamics (S h -I h -S m -I m ) 


100 


50 100 


150 200 250 300 350 

Day 


2 

0 l-------■— 

0 50 100 150 200 250 300 350 

Day 



15 


10 


0 50 100 150 200 250 300 350 

Day 


Figure 5: Population dynamics of the honeybee-mite-virus model (7) when r = 1500; K = 160001; p = 0.9; <4 = .15;/x^, = 
0.1; a = 0.05; dm = 0.1;/im = 0.01; c = 0.005; /3 h = .3 ■ p mh = 0.08; p mh = .001; $ hm = 0.03 and S h (0) = 410,4(0) = 
10, 5 m (0) = 35, 7 m (0) = 10. The healthy honeybee population Sh is in red; the virus infected honeybee population 4 is in 
black; the healthy mite population S m is in blue; and the virus infected mite population 7 m is in cyan. 


model by using the susceptible-infectious (SI) disease framework. Our proposed model includes the parasitic 
interaction between Varroa mites and honeybees as well as different virus transmission modes occurring 
at different phases of Varroa mites. More specifically, we use frequency-dependent transmission functions 
to model horizontal virus transmissions between Varroa mites and honeybees during the phoretic phase of 
mites and use a Holling Type I functional response to model parasitic interactions between Varroa mites 
and honeybees. We also use a density-dependent transmission function to model the vertical virus trans¬ 
mission between Varroa mites and honeybees during the reproductive phase of Varroa mites. We also apply 
the frequency-dependent transmission function to model the horizontal virus transmission among honeybees. 

Summarize the main findings: Our analytical and numerical results of the full system suggest the 
following: 


1. Initial honeybee populations play an important role in its persistence since its dynamics exhibits strong 
Allee effect in the absence of both parasites and virus (see Theorem 3.1). In addition, patterns of popu¬ 
lation dynamics are sensitive to initial conditions as suggested by our numerical results, i.e., depending 
on initial conditions, the full system can experience the catastrophic event where honeybees collapse 
dramatically, or both mites and honeybees can coexist but exhibit different dynamical patterns. 

2. In the absence of Varroa mites, the honeybee and disease system has only equilibrium dynamics (see 
Theorem 3.4 and the first row of Figure 2-3). In the presence of Varroa mites, the synergistic effects 
make disease and mites more persistent (see Theorem 3.5-Theorem 3.6), and thus, difficult to control 
(see Theorem 3.2 and the related arguments). In addition, the synergistic effects make disease and 
mites can also drive both honeybees and mites go extinct (see Figure 5). 
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3. In the absence of virus infections, the Varroa mite and honeybee system can be destabilized by the 
low adult workers to brood ratio in the colony where the system has oscillating dynamics leading to a 
sudden extinction of all species (see Theorem 3.3 and Figure 1). This dynamical property is called the 
catastrophic event which has been inherited by the full system when the virus infection is present (see 
Figure 5). This phenomenon could be linked to CCD which has been observed in honeybee colonies 
where low hive to brood ratio could be a contributing factor. 

4. Our numerical simulations suggest that large values of virus transmission rates can drive the extinction 
of healthy mites while extremely large values can lead to the extinction of all species. 


Identify the contribution to the broader field: Our theoretical results in Theorem 3.1 imply that 
the persistence of honeybees requires a critical threshold of the adult workers to brood ratio. This illustrates 
previous observations that a proper adult worker to brood ratio is needed to have a successful honeybee 
colony [46] . Our current work indeed provides many useful insights on our proposed questions regarding the 
complicated synergistic effects of Varroa mites and the associated virus infections on the honeybee popula¬ 
tion dynamics and the persistence as well as how these synergistic effects may potentially contribute to CCD 
of honeybee colonies. Need more help from Gloria 

Compare the findings with other works: As we mentioned in our introduction, the work of Sumpter 
and Martin [54] and Ratti et al. [42] are most relevant modeling papers for our study proposes. However, 
Sumpter and Martin [54] assumed that the mites’ population is constant and the virus transmission occurs 
only through Varroa mites. This assumption prevents us to study how mites’ population affect the virus 
transmission and the honeybee population dynamics. The model of Ratti et al. [42] does not include the 
virus transmission through honeybees themselves which leads to the extinction of virus when Varroa mites 
go extinction. This is not realistic. Our honeybee-mite-virus model allows us to investigate the honeybee 
population dynamics in the absence of both the Varroa mites and virus (the honeybee dynamics (11)), in the 
absence of the virus (the virus-free dynamics (10)), in the absence of the Vorroa mites (the mite-free dynam¬ 
ics (8) , where the virus is transmitted through other contacts), and in the presence of both Varroa mites and 
virus to explore the synergistic effects of Varroa mites and the associated virus infections on the honeybee 
population dynamics as well as the potential mechanisms contributing to CCD of honeybee colonies. More 
specifically, our model is able to provide sufficient conditions that can guarantee the persistence of honeybee 
population and conditions that lead to the extinction of honeybees. 

We would like to mention the recent work by Perry et al. [41] which examined the social dynamics 
underlying the dramatic colony failure with an aid of a honeybee population model. Their model does not 
include population dynamics of Vorroa mites or virus infections, but it does includes bee foraging perfor¬ 
mance varying with age, and displays dynamics of colony population collapse that are similar to field reports 
of CCD. We plan to expand our current modeling framework to include the detailed stage-structure model 
of honeybees with social dynamics. 

Even if nutrients are not limiting, as assumed in the VARROAPOP model, colonies go to extinction due to 
reduced longevity of adult workers, adult population decline, and a reduced ratio of adults to brood. This 
ultimately reduces brood rearing so that bees die at rates that are greater than emergence of new adults. 
This begins the downward spiral in the adult population that ultimately leads to colony death. 

Suggest options for further research: In our proposed model, we assume that a fixed ratio of adult 
workers to brood which does not reflect the reality since this ratio varies with the availability of nutrient 
that depends on the quality and availability of pollen. The mortality of honey bees is resulted from differ¬ 
ent risk factors such as parasites, pathogen, viruses, pesticides, nutrition and environmental changes [52]. 
Among these factors, honey bee nutrition is one of the most fundamental factors which impacts honey bee 
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health and influences the capabilities of honey bees to combat different stressors. Nutrients such as pro¬ 
tein from pollen is essential in fighting parasites and diseases, and in maintaining the high adult worker to 
brood ratio in honeybees. This is because that honeybees solely depend on pollen and honey for satisfying 
their daily needs where pollen provides protein, essential amino acids, vitamins and mineral for bees while 
honey or nectar are their carbohydrate resources [3] . Cage studies showed that honey bees could survive for a 
long time without pollen provision [5], but pollen feeding significantly prolonged their lifespan [29, 34, 49, 48]. 

In addition, intensive research has investigated the connection between nutrition and honey bee disease 
and stress resistance. Both cage studies and field studies indicated that bees with poor nutrition were under 
more stress (Wang et al. unpublished data) [60], more susceptible to Noserna and Varroa destructor, and 
had shorter lifespan [17, 43]. DeGrandi-Hoffman found a similar result that bees with good pollen or protein 
nutrition were more resistant to many different viruses [13] . Furthermore, studies on molecular mechanisms 
suggested that pollen nutrition may positively affect antimicrobial peptides and improve immune defensive 
response to parasites [15]. Therefore, a more realistic model that includes nutrient/brood dynamics, explicit 
division of labor, and seasonal effects, are needed. This is our on-going project. 


5. Proof 

Proof of Theorem 3.1 

Proof. According to Theorem A.4 (p.423) of Thieme (2003), we can conclude that Model (7) is positive 
invariant in A'. Let d = mm{dh , d m }, Nh = Sh + Ih, N m = S m + / m , and N = cNh + N m . Then we have 

N' = cN' h +N' m = ~ cd h N h - cnhl h - d m N m - » m I m 

< K-lwc) 2 ~ min { dh ’ d m}(cN h + N m ) = - min{4, d m }N = J^ N2 - dN 


r_|_* / fr.\ 2 —4 

which implies that lim sup^^ N(t) < c— — v ^ - = N* with implication that lim sup^^ N(t) = 0 if 

either ^ < 2 \/~K or IV(0) < N c = holds. The arguments above also imply that if ^ > 2 \[k, 

then we have 


lim sup Nh(t) < N*/c and lim sup N m (t) < N*. 

t—> oo t—>oo 

Similarly, we have follows for Nh'. 

K=S' h +I' h = - d h N h - Hhh - aN h N„ 


< 


K+(S h +pI h ) 

rN 2 h 
K+Nl 


— dhN h 


which implies that 


+ 


lim sup Nh (t) < 

t—} OO 


(*)' 


-AK 


= N* h if—> 
dh 
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and 


lim sup Nh (t) = 0 if -j- < 2 \f~k or Nh( 0) < Sf = 
dh 
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On the other hand, we have 

N'h = S' h + I' h > ~ (dh + Ph + aN m )N h 


> $$s-(d h + iM h + aN*)N h 


> 


r Nf 


K/p 2 +N 2 

Therefore, apply the comparison theorem, we can conclude that 


— {dh + Ph + cxN*)Nh 


lim inf Nh (t) > 

t—¥ OO 

if the following inequalities hold 


dh+Ph+aN* 


{( 


dh+Ph+aN * 


- 4A7p 2 


=m 


dh + Ph + cxN* 


K d h +ii h +aN* 

>2 \l— and N h (0)>NZ = 
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The discussions above provide sufficient conditions that allow Nh being persistent, i.e., when the inequalities 
d h +h h +aN• > 2-y/p- and Nh{ 0) > N_h hold. This implies that Sh is also persistent under this condition since 
all species go extinct if Sh goes extinct. However, this persistence condition does not provide an estimate of 
Sh- To explore an estimation of Sh, we look at the population of Sh- 


Recall that if dh+fi y +aN , >2y and Nh{ 0) > N%., then we have the following inequalities 

N * < lim inf Nh{t) < lim inf N{t)/c < lim sup N{t)/c < N*/c. 
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which implies that 


dh+Ph+ ^y* +0mh+oc)(N*-cNl) \ y d h +0 h + +(Pmh+<x)(N*-cN* h ) 


-4 K 


lim inf Sh ( t ) > St = 

t—> OO 

Therefore, we can conclude that Sh is persistent with the following properties: 

< lim inf Sh{t) < lim inf Nh(t) < lim inf —— < lim sup N(t)/c < N*/c 

t—t oo t—> OO t—> OO C t > CO 


if the following inequalities hold 

1. > max \d h +p h + + 0 mh + a)(N* - N£), d ^+ aN ' } with JV h (0) > S h ( 0) > S c h . 

It is easy to check that the extinction equilibrium E 0 = (0,0,0, 0) is always local stable. We omit the details. 
Based on the discussion of the upper bound of the total population N and the total population of honeybees 
IV/,, we can conclude that the system (7) converges to E 0 globally if 2 \flc > -ff- holds; while if the initial 
population satisfies either IV(0) < N c or IV/, (0) < Nf, then the system (7) converges to Eo locally. 

□ 


Proof of Theorem 3.2 

Proof. Now we consider the population of N m . From Model (7), we obtain the following inequalities: 

cotNhNm ( d m ~t~ /r m )IV m ^ N m — S m + I rn — caNhNm d m N m I rn IV m (co Af, dm) 

which implies that if N* < then we have 

N'm < N m [caN h - d m \ = N m [ca{N - N m )/c - d m ] 

< N m [aN* -d m - aN m } < N m [aN* - d m ] 


This indicates that limsup^^ N m (t) = 0. 


Assume that the following inequalities hold > 2 \J~K, 
3.1, we have 


and Sh{ 0) > S^. Then according to Theorem 


lim sup Sh (t) < lim sup Nh (t) < N£. 

t—t OO t—y OO 


Now let us focus on the population of S m . Notice that we have the following inequalities when time is large 
enough, 


q> 


= S n 


q _ Phmll i _ J 

eas h Sh+Ih a„ 


— Sm [cCViV/, dff! ' Sm \caN h dm\ 


This implies that if N% < ^ and -ff- > 2 \[k, and S/,(0) > Sf, the healthy mite population S m goes extinct 
while the total population of honeybees persists, i.e., 


lim sup S m (t ) = 0. 

t—> OO 

Now we look at the population dynamics of I m and //,. Let / = clh +/ m . From Model (7), then we have 
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the following equations: 
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Therefore, if 
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> j^min ^/3 h ,cj3mh + cfimh + caSl } - max j (dh + Mh), {d m + Pm )}] > 0 

which implies that the disease I persists by applying the average Lynapunov theorem (Hutson 1984). 
On the other hand, we have the following inequalities: 
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Pp- + cpmh + caSh } - min j(t4 + fi h ), {dm + Mm)}] 
+ cfimh + caSh } - min | {dh + Hh), {dm + Mm)}] 


Assume that 
inequalities: 


> 


2 \J K 


d h +hh+aN , N h {0) > 7V£. Then according to Theorem 3.1, we have the following 


N t < lim inf Nh {t) < lim sup Nh (t) < < lim sup N{t)/c= N*/c. 


£—>■00 


£—>■00 


This implies the following inequality when 


max \Ph+ ^ h ™y iCPmh+cPmh+caN^ 


ain ■} 


min ( dh+hh),(dm+hm) 


} 


i < 1, 


then we have 


iL 

1 


< max 


\i h =i m =o 

which implies that the disease goes extinct, i.e., 


jdh + . C Pmh + cfimh + C0tN£ } - min | {dh + Mh). {dm + Mm)}] < 0 


lim I{t ) = 0. 
£—>00 


□ 
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Proof of Theorem 3.3 

Proof. If M = 0, the virus-free subsystem (10) reduces to the only healthy honeybee population: 


Si = 


r q2 

h - d h S h 


K + Sl 

which leads to the following three boundary equilibria if 

2K > dh: 

(0,0), (N£, 0), and (IV,*, 0). 

The local stability can be easily determined by the eigenvalues evaluated at its Jacobian matrix. Simple 
algebraic calculations show that (iV£, 0) is a saddle if Nf < while it is a source if Nf > and (N£, 0) 
is a sink if NT < — while it is a saddle if NT > 4m.. 

tl OLC <1 OiC 

Now let (H,M) be an interior equilibrium of the virus-free subsystem (10), then we have the following 
equations hold: 

0 = caHM - d m M O H = 4m 

r,L ca 


o = 


= - “™ « M = s - 4 ) 

which gives the unique interior equilibrium (H *, M *) = ^ ^ 


rH * 


k K+(H *) 1 2 


- d h )) provided that Nf < 4m < 


N* h - 


The local stability of ( H*,M *) is determined by the eigenvalues A i,i = 1,2 of the following Jacobian 
matrix of (10): 


Jhm — 


( rH*(K-(H*) 2 ) 
(k+(H*) 2 Y 


- olH* 


\ 


acM* 


which gives the following two equations: 

Ai + A 2 


A 1 A 2 


r H*(K-(H*) 2 ) 

{k+(H-yf 

ca 2 H*M* 


( 12 ) 


(13) 


Therefore, we can conclude that ( H*,M *) is a sink if H* > \/~fc while ( H*,M *) is a source if H* < \Z~K. 
The discussion above implies follows: 

1. If —'-= < dh , then the extinction equilibrium (0,0) is the only locally stable equilibrium for the 

2 V k 

virus-free subsystem (10). Thus it is globally stable. 

2. If ^j— > dh and either IV f or ^ > Nf. then the virus-free subsystem (10) has three boundary 

equilibria: (0,0), (5£,0), (N£,0) where (0,0) is always locally stable. 

• If Nf > 4m, then (S^. 0) is a source and (1V^,0) is a saddle. Since (10) is a two-D ode system, 
it has global stability at the extinction equilibrium (0,0) according to the Poincare-Bendison 
theorem [40]. 
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• If N* < then (S£,0) is a saddle and (N£, 0) is a sink. This implies that (10) has two locally 
asymptotically stable boundary equilibria (0,0) and (N%, 0) which are reserved as the only two 
attractors for the model. 

3. If > dh and N% < < N%, then (10) has three boundary equilibria (0,0), (S£,0), (N£, 0) and 

the unique interior equilibrium ( H*,M*) where (0,0) is locally stable, both (S£,0) and (N£, 0) are 
saddle nodes. The local stability of (. H*,M *) is determined by the sign of H* — \[K: if H* — \[k > 0, 
then ( H*,M*) is locally asymptotically stable while if H* — \f~fc < 0, then ( H*,M *) is a source. 

Define Rq 1 = . Then we can conclude that (10) has no interior equilibrium if Rq 1 < 1 or *o M > & > l; 

and (10) has a unique interior equilibrium (H*,M*) if 1 < Rff < pf. More specifically, ( H*,M*) is locally 
stable if 1 < Rff < Sh- while it is a source if Rff > max |l, ^□ 


Proof of Theorem 3.f 

Proof. If / = 0, the mite-free subsystem (8) reduces to the only healthy honeybee population: 


Si, = 


r c2 

h - d h S h 


K + Sl 


which leads to the following three boundary equilibria if -pp > dh- 


(0,0), (N£, 0), and (IV*, 0). 


The local stability can be easily determined by the eigenvalues evaluated at its Jacobian matrix. Define 
Rq = ■ The simple algebraic calculations show that if IRj 7 < 1, (Nfc, 0) is a saddle and (N£,0) is 

a sink; while if Rq > 1, then (IV)), 0) is a source and (N%, 0) is a saddle. According to Theorem 3.1, the 
extinction equilibrium (0, 0) is always a sink. 


Let (S, I) be an interior equilibrium of the mite-free subsystem (8), then it satisfies the following equations: 


— (dh + Ph) — 0 


r(S+pI ) 2 , o __ q Phi _ n 
K+( S +pip ah * * S+I - U 


S = , ^ J = —J- -= a l 


Ph~(dh+Ph) 


W 

d h,+! J -h 


rl 


_ K , ; J2 

(p+p ) 2 


(dh + 


-l 


= 0 


r{S+pI ) 2 _ j q _ q Phi _ n 
K+(S+pI ) 2 h *s+/- u 


^S±tnL-d k s-(d h + n )i = o 


(14) 


rl 


-—— - ’T -\~I 2 

(» + P)‘ 


— ((a + 1 )dh + Ph) — 


rl 


—^r+/ 2 

(p+p > 2 


-d = 0 


T+a = d h + hh- 


Therefore, if Rq = > 1 and £ > then we can conclude that the mite-free subsystem (8) can 

have two interior equilibria (S^,/^),fc = 1,2 where 


1 


Ph 


P h 




~ 1 


d = (a + 1 )d h + Ph = d h 


dh+Ph 


Ph 


. d h -\-p h 


- 1 


Ph, 


30 
























and 


Ih = 




A K r 

^ 4 ++tf , d 

—, I h = - 


(5)' 


/I K 

4 f+tf 


, = alfr, k = 1,2. 


2 n 2 
Now we examine the local stability of these two interior equilibria provided that they exist. The Jacobian 
matrix of (8) evaluated at the interior equilibrium ( S , I) = ( al , I) can be expressed as follows: 


/ 2rK(a+p)I _ , _ 0 h 2rpK(a+p)I _ a 2 ^ fe 

((a+p) 2 I 2 +K) 2 h (1+a) 2 ((a+p) 2 / 2 + A') 2 (1+a) 2 


Jhv — 

whose two eigenvalues A,;, i = 1,2 satisfy the following equations: 


V 


0h 

IT+^F 


a0h 

■(T+^F 


J 


Ar + A 2 = „ 2rK tt p) L -d h --^~ 


2K 


r(a+p) 2 I 


((a+p) 2 I 2 +K) 2 


_ 2K _ 

(( a+p) 2 I 2 +K)(a+p)' 


1+a {{a+p) 2 1 2 +K) (a+p) ( a+p) 2 I 2 +K 

(dh + - d h - rh = {d h + 


-d h --& 


1+a 


2aK 


_((a+p) 2 I 2 + K)(a+p) 


- 1 


< 




2 K 


(q+p) 2 I 2 +K 


- 1 


and 


\ \ _ Ph 

a i a 2 - (T+^F 


0h 


+ adh 


q0h 


(1+a) 2 


2 raK(a+p)I 
' ((a+p) 2 J 2 +K) 2 

2rA~(a+p) 2 J . , a/+ 

((a+p) 2 / 2 +A) 2 ^ UUh ^ 1+a 


2rp.K’(a+p)I 


q 2 0h 


(1+a) 2 ((a+p) 2 / 2 + A) 2 W (l +a )^ 


1 Ph 

2 rk(a+p) 2 I . 7 

\ - (1+a) 2 

((a+p) 2 / 2 +lC) 2 U J 


If J Iq = > 1 and ^ then we have two interior equilibria (S*, +), fc = 1, 2 where 


d a+p 


4 = 




_ A K r 

4 (a+p) 2 r2 3 


( 5 )' 


j A- 
4 (a+p) 2 


2 r ‘ 2 
This implies that Ai(+)A 2 <+) < 0, Ai(i|)A 2 (/ 2 ) > 0 and 


, Su = alu,k = 1 , 2 . 


Ai(^) + A 2 (/2) < 


< 


(dfc 


2K 


Ph 

1+a 


_(a+p) 2 I 2 + K 


2K 


— 1 


- 1 


'++p) 2 ~ 


= 0 


(15) 


(16) 


(17) 


since 




r 


> —= > 


2d 


Vk 

(a + p) 


=► U") 2 > 


k 

+++■ 


Therefore, we can conclude that if (Rq = > 1 and ~ , then the system has two interior equi¬ 

libria (S%, +), fc = 1,2 where (+, ID is always a saddle while (S' 2 ,1 2 ) is always locally asymptotically stable. 


Notice that ~ > 

a 


2 \fk ,. i- 
a+p 2\fk 


> 


d 

q+P 


> dh- The discussions above show that the 
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mite-free subsystem (8) has no interior equilibrium if either 3?^ < 1 or holds. This leads to the 


following two cases: 
1. If either 


2v K 


< dh , then the mite-free subsystem (8) has only the extinction equilibrium (0,0) which 


is locally stable, thus it is globally stable. 


2. If > 1 and d h < 
equilibria 


2 yJK 


< 


a+p 


then the mite-free subsystem (8) has the following three boundary 


(0,0), (N£, 0), and (N £, 0) 

where (0, 0) is locally asymptotically stable; {Nf, 0) is a source; and (N£,0) is a saddle. Then according 
to the Poincare-Bendison theorem, we can conclude that all interior points of converges to the 
extinction equilibrium (0,0). 


3. If IRq - " < 1 and dh < 


2V K 


then the mite-free subsystem (8) has the following three boundary equilibria 


(0,0), (Nl 0), and (iV^,0) 

where (0,0) is locally asymptotically stable; (Nfc, 0) is a saddle; and (iV^,0) is a sink. This implies 
that both (0,0), and (N£, 0) are locally asymptotically stable. 


□ 


Proof of Theorem 3.5 

Proof. Let Nh = Sh + Ih, then from Model (9), we obtain 


K 


r(Sh+p!h ) 2 

P~\~(Sh~^pIh) 2 


d h N h l^hlh 


< 


k+Ni 


— dhNh 


which implies that when 


2 yfk 


> dh, we have 


limsup (Sh(t) + I hit)) = limsup N h (t) < 

t—yoo t —^oo 


r 

dh 



= N* h . 


Therefore, we can conclude that if the inequalities ^ > 2 \J~k. and < d m + [i m hold, then we have 
limsup^^ S m (t) = 0 since 


IL 


I | g {dm, + Mm) 

— Ca^ 

Y„ {dm + Mm) 

< calm 

Y* {dm + Mm) 

COL 


ca 


ca 


This implies that the limiting dynamics of Model (9) is reduced to the mite-free model (8). 


An interior equilibrium (Sh, Ih, Im) of the healthy-mite-free subsystem (9), satisfies the following equa- 
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tions: 


COilrr 


Ih + s h ~ ^4^ 


= 0 =>N h = I h + S h = **±± 


- AJV * ~ "*'■ - » 7 * 


0 ^ k+TsT+p^F dhNh aNhIm P hIh ~ 0 

kXirXZi - dM ‘ ~ aN ^ - “ A = 0 


\+p i hf 


Im. — 


*+( 


dm ~\~^n 


h+pi> i) 




( dm+Hm j 


= fi{h) 


(18) 


S, 


Phlh | Pmhlm | O T 
Su + Tu ' Su+Tu ^ 


S h +I h T Sh +I h 

(^-4) - 


dm+flrr 


CCtfimhlrr 

dm+4m 


fimhlrr 


&lhlm H“ (dh H“ H , h)dh 
(xdhlm H“ (dh + fd'h'jlh 


I'm. — 


Im. — 


( Sa ^ BL - Ih ) 


,+ J m +i +Pmh\Ih- 


. dm+fj-n 




= /a(4) 


The equations above imply that the interior equilibrium of (9) is the positive intercept of fi(Ih ) and f 2 {Ih) 
subject to 0 < Ih < m '4' m • The expression for the function fi(Ih) implies that the subsystem (9) has no 
interior equilibrium if 

dh{dm + Mm) j / j . \ 

■v 1 ^ TC(X ^ dh\d m + Mm) 


r < 


since / m = fi(Ih) < 0 when this inequality holds. 


Assume that 


2v K 


> rfh +^".+aA and ./V^O) > IVjj, then according to Theorem 3.1 and 3.2, then we have 
N!h < MX — limiaf Nh(t) < limsup Nh{t) < 


which implies that the set {Sh + Ih > N_h} is invariant. If I m = 0, then the subsystem (9) reduces to the 
mite-free system (8). According to Theorem 3.4, the omega limit set of the mite-free system (8) is (N£,0) 


when 5Jn = 


Ph 


0 d h +n h 


< 1 holds. If Ni > dm+tlrn then we have 

rl r.rv 1 


felr m =o = c «(^-^)> 0 ' 


This implies that the disease I m persists by applying the average Lynapunov theorem (Hutson 1984). Notice 
that 

T I _ PmhS h I m > Q 

h \l h =0 N h > U - 

Therefore, the virus infected honeybees Ih also persists. 

□ 


.(19) 
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Proof of Theorem 3.6 

Proof. The first part of Theorem 3.6 can be deduced directly from Theorem 3.1, 3.2, 3.3, 3.4, and 3.5. We 
focus on the sufficient conditions that lead to no interior equilibrium (Item 3) and the persistence of mites 
(Item 4). 


If (Sh, Ihi Smi Im) is an interior equilibrium of the system (7), then it satisfies the following equations: 
0 = S, 

caSh—d■ 


O _ Phmlh _ J 

cas h Sh+Ih a„ 


=> h=S h 


dm+Phm—caSh ’ ca 


< s h < 


( 20 ) 


^ = S H + I H = ^ + ^ ea s h =g 1 {S H ) 


0 = ^^^s-dh{S h + I h )-a(S h + I h )(S m + I m )-ti h Ih 


-g^^j-d h lS h +I K )-^ h I h 
AT _ C ,T _ K +( s h +Plh) _ 

JV m - Om+lm - a(S h +I h ) 


( 21 ) 


jy _ rS h (0 hm + (d m -caS h ))(j3 hm + (d m -caS h )(l-p)) 2 _ dh _ Ph(caS h -d m ) _ (Q \ 

m aiphm \ K{hhm+(d m — caSh)) 2 + St l (idhm + (d m -caSh)(l-p)) 2 ^ “ Phm 


0 — COc(Sh T hi.) {Sm T Im) (I m ( S r n. T Im) /dmln 


j N rn (ca(Sh-\-Ih) — dm.) \cafihrn.Sh — dm (hhm + (^m — COlSh ) ) ] N m (caSh—d m ) (Phmp-dm) 

V™. Pm($hm+(dm,-caS h )) p. m (Phm + (d m -CaS h )) 


( 22 ) 


= N m g 3 (S h ) = g 2 (Sh)g3(Sh) 


0 — Sh ) h + "b Anh-lmJ — OcIh(Sm + Im) ~ (dh + gh)Ih 


dm+Phrn~ caS h 

caSh — d m 


Ph( ca S h -d m ) . o N m (caS h -d m ){P hm +d m ) 
r omh 

Phm PmPhmSh 


I 5 Nm(caS h -d Tn )(f3 hrn +d m ) 
m pm(Phm + (d m -caS h )) 


- aN m - (d h + gh) = < 


Ph( d m+Phm-c<xS h ) _j_ P mh N m (P hm +d m )(d m + P hm -caS h ) P mh N m (P hm +d m ) _ „ \ _ 

Phm. PmPhmSh d-m \ H 1 ) 


Nm = 


d h (dm + P h m-O^S h) 

Phm 


. Pmh(Phm~ >rd m)(.d rn +P hm _ ca.S h ) _ 0 rnh (0 hrn -\-d rn ) 

PmPhmSh Prn 


= gi(Sh) 


• (23) 
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Therefore, the interior equilibrium (Sh, Ih, S m , I m ) are positive solutions of the following four equations: 


d m +p hm -c aSh - " ca 

[c a s h -d m )0 hm +d m ) < 1 subject t0 4% < s h < dm+ J hm , 

Mm — caS^ J COL COL rfm+fem+Mm 

_ rSh(Phm + (dm-coiSh))(Phm + (dm-co:Sh)(l-p)) 2 _ d h hh(. ca S h -d m ) 

a ~P hm [k(p hm +(d m -c a S h )) 2 + Sl{p hm + (d m -cccS h )(l-p )) 2 ] “ Phm 

P h (d m +9 hm -c a S h ) + 

_ Phm _ 


N h 

— Sh ~\~ dh — 

Im 

= gs(Sh) = - 

N m 

N m 

= g2(s h ) = - 
< 

N m 

= g*(S h ) = - 


(24) 


Pm.h(Phrn.+ d rn)(drn+Ph.m.- co ‘ S h) Pm.h.(Phm.+ d m) 
Pm0hm S h Aim 

MmSfc [Ph( d m+Phrn)-Phm( d h+Vh)- coc Ph S h\ 


[cxVmPhrn+ cc *Pmh(P h m+dm)-PmhPh r n(Phm+dm)]Sh-Pmh(Phm+dm)(d rn +/3hm) 

Thus, we can conclude the following statements regarding the sign of <74 (S^): 

1 . If Phidm T Phm) ^ Phm(dh T 9h) and QLjJjmPhm T CQ!.P rn h(Phm T dm) ^ Pmh Phm (Phm T d m ), then 
94 (Sh) > 0 for all S h > 0. 

2- If Ph{d-m T Phm) ^ Phmidh “1“ 9h) ^-Ild OtHmPhm T CO.P 7 nh{Phm “1“ ^m) ^ Pmh Phm {Phm T d m ), then 
94 (Sh) > 0 when 


Pmh{Ph m + d 

m )(d 

m + Phn 


a[l m Phm T CQ^Pmh^Phm T d rn ) PmhPhm(Phm T d n 


< s h < 


Ph(d m + ph m ) ~ Ph m (dh + Hh) 
caph 


3- If Ph(dm T Phm ) < Phm(dh ~h [dh) and CX[lmPhm T Caf3mh(Phm ' dm) > PmhPhm(Phm T dr n ). then 
gp(Sh) > 0 when 

Pmh (Ph m + d 

m )(d 

m + Ph m) 


s h < 


& 9 m Phm T C(xPmh(.Phm T d m ) PmhPhm(Phm T dn 


4. If Phipm T Phm) Phmidh T and OiflmPhm T CC^Pmh i^Phm T d m ) <C Pmh Phm (Phm T d m ), then 
^ 4 ( 5 ^) > 0 when 

m + Ph m )-Ph m (dh + [dh) 


S h > 


caph 


Notice that the interior equilibrium (S h , I h , S m , I m ) requires g 3 (S h ) > 0 , i.e., ^ < S h < dm ^ hm - 
Therefore, the interior equilibrium (Sh, Ih, S m , I m ) does not exist if one of the following inequalities hold 

1 . > - Af-, > p m J hm - ca and 

Ph _ Phm 


2 . 


3. 


dh^'h'h 

> 

dm+Phm 

7 Phm~\-dm 




Pmh 


CX-hm (Phm 

T caPmh 


(d 

m+fihm ) 

p h 

< 

Phm 

ap m Phm 

dh+Ph 

dm ~\~ fihm 

’ $hm~\~dm 

Ph 

> 

Phm 

aPmPhm 


> 


*m ~r~ 9m 


^ Pmh Phm and 


or 


$mh(dm+fi h r, 


dh+Vh dm+Phm / j 

- ^ tlm • 


a [dm Phm 
( d m+Phm ) 


+cap mh -p mh p hn 


< dm 
ca 


Ph. Phm. 

_ V (i , fl, _ , YV arlr l d h+^h dm+Phm {drn+JMn ){dm+Phn 

d »+ldh ^'-WmiphS' PmhPhm 0h(d m +0 hm )(d h +tt h ) d m +Phm+Pm 
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Assume that —> dh+llh + aN and Nh( 0) > N)., then according to Theorem 3.1 and 3.2, then we have 


N-h < NX — liminf Nh(t) < limsup Nh(t) < 

i-^ 00 t-y oo 


which implies that the set {Sh + Ih > NX} is invariant. If S m = I m = 0, then the full system (7) reduces 
to the mite-free system (8). According to Theorem 3.4, the omega limit set of the mite-free system (8) is 
(N£, 0) when 01q = dh ^ llh < 1 holds. If > dm +^ m , then we have 

fXmJ > ca(N? - dm+AIm ^ > 0 

N m \ N m =0 ~ Ca \ h ca ) ^ U ' 

This implies that the mite population N m persists by applying the average Lynapunov theorem (Hutson 
1984). 
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